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Abstract 

I analyse a generalised Burger's equation to develop an accurate 
finite difference approximation to its dynamics. The analysis is based 
■ upon centre manifold theory so we are assured that the finite difference 

Q\ . model accurately models the dynamics and may be constructed sys- 

tematically. The trick to the application of centre manifold theory is 
£>V to divide the physical domain into small elements by introducing insu- 

i ■ lating internal boundaries which are later removed. Burger's equation 

is used as an example to show how the concepts work in practise. The 
^ (— i resulting finite difference models are shown to be significantly more 

O ■ accurate than conventional discretisations, particularly for highly non- 

linear dynamics. This centre manifold approach treats the dynamical 
k>( ; equations as a whole, not just as the sum of separate terms — it is holis- 

$_i ■ tic. The techniques developed here may be used to accurately model 



the nonlinear evolution of quite general spatio-temporal dynamical 
systems. 
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1 Introduction 

I introduce a new approach to finite difference approximation by illustrating 
the concepts and analysis on the definite example of a generalised Burger's 
equation. In some non-dimensional form we take the following partial differ- 
ential equation (pde) to govern the evolution of u(x,t): 

du ^ du d 2 u 3 
dt dx dx 2 



XXI 



This example equation, which includes the mechanisms of dissipation, u 
and nonlinear advection/steepening, auu x , generalises Burger's equation by 
also including a nonlinear damping, j3u 3 . Consider implementing the method 
of lines by discretising in x and integrating in time as a set of ordinary 
differential equations. A finite difference approximation to ([]]) on a regular 
grid in x is straightforward, say Xj = jh for some grid spacing h. For example, 
the linear term 

d 2 u _ Uj+x - 2 Uj + uj-i Q m 



dx 2 h 2 

However, there are differing valid alternatives for the nonlinear term uu x : 
two possibilities are 



du _ Uj(u j+1 - uj-i) / 2 n _ u] +1 - u 2 



u 



dx 2h V ) 4h 



Which is better? The answer depends upon how the discretisation of the 
nonlinearity interacts with the dynamics of other terms. The traditional 
approach of considering the discretisation of each term separately does not 
tell us. Instead, in order to find the best discretisation we have to consider 
the dynamics of all terms in the equation in a holistic approach. 

Centre manifold theory (see the book by Carr |IJ) has appropriate char- 
acteristics to do this. It addresses the evolution of a dynamical system in 
a neighbourhood of a marginally stable fixed point; based upon the linear 
dynamics the theory guarantees that an accurate low-dimensional descrip- 
tion of the nonlinear dynamics may be deduced. The theory is a powerful 
tool for the modelling of complex dynamical systems [0, H, O, O, 0, e.g.] 
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such as dispersion |23], [TI], |26], e.g.], thin fluid films ||, [T7|, [H], e.g.] and other 



applications discussed in the review | 19fl . Here we place the discretisation of 



a nonlinear pde such as (|1]) within the purview of centre manifold theory by 
the following artifice (such mathematical trickery has proven effective in thin 
fluid flows Hl7|). Introduce a parameter 7, < 7 < 1: at the midpoints of 
the grid, x = (j + l/2)h, insert artificial boundaries which are "insulating" 
when 7 = 0, but when 7 = 1 the boundaries ensure sufficient continuity to 
recover the original problem. In essence this divides the domain into equi- 
sized elements centred upon each grid point, say the domain is partitioned 
into m elements. For (|Tp we may use 

du + du~ , , h ( du + du~ \ / \ 

where u + is just to the right of a midpoint and u~ to the left. When 7 = 1 
these reduce to conditions ensuring appropriate continuity between adjacent 
elements. When 7 = they reduce to conditions equivalent to the insulating 

du + du~ 
dx dx 

We treat terms multiplied by 7 as "nonlinear" perturbations to the insulated 
dynamics. Then in the "linear" dynamics governed by 

du d 2 u 



dt dx 2 ' 

each element evolves exponentially quickly (in a time O (h 2 )) to a constant 
value, say u = Uj in the jth element (the particular value depends upon the 
initial conditions). But in the presence of the perturbative influences of the 
nonlinear terms and the coupling between elements, the values Uj associated 
with each element will evolve in time. This picture is the basis of centre 
manifold theory which is applied in Section ^| to assure three things: 

• the existence of an m dimensional centre manifold parameterized by 

Uj] 

• the relevance of the m dimensional dynamics as an accurate and sta- 
ble model of the original dynamics (0) (sometimes called asymptotic 
completeness [p2| or 0); 



and that we may approximate the shape of the centre manifold and the 
evolution thereon by approximately solving an associated pde. 
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These dynamics on the centre manifold form a finite difference approxima- 
tion. For example, the analysis in Section |3] of the generalised Burger's 
equation (|1|) is unequivocally in favour of the discretisation 



iij + auj 



2h ~V2 V ' 



- 2uj + Uj-t) 



Uj + i — 2Uj + 



h 2 

- 7^ (wf-i + Suj-iUj + 16u 3 j + 3M j+ iM| + , (3) 

as an early approximation (iij denotes duj/dt). Provided the initial con- 
ditions are not too extreme, centre manifold theory assures us that such 
a discretisation models the dynamics of ([I]) to errors (9(||m|| 4 ). Observe 
that it is best to discretise the nonlinear uu x terms directly, but that there 
is a nonlinear correction involving the second difference. Further, the cu- 
bic nonlinearity is discretised in a non-obvious complicated form. Margolin 



& Jones [10 1 have previously applied inertial manifold ideas to discretisa- 
tions of Burger's equation. However, they used just two basis functions on 
each element and invoked the adiabatic approximation for time derivatives 
of "slaved" modes. Here I allow arbitrarily convoluted dependence within 
each element and include all effects of time variations. The methodology 
presented here provides a rigorous approach to finite difference models. 

The discretisation (|3p is just a low-order approximation, centre mani- 
fold theory provides systematic corrections. For example, Equation (|3]) is 
obtained from terms linear in the coupling 7. Analysis including quadratic 
terms in 7 leads to the fourth-order accurate model discussed in Section |j. 
Analysis to higher orders in the nonlinearity, discussed in Section ^ shows 
higher order corrections to the discretisation of the nonlinear terms and also 
incorporates effects from a coupling of the different nonlinearities. 

Computer algebra is an effective tool for modelling because of the system- 
atic nature of centre manifold theory |I4| , |8|, p0| . The specific finite difference 
models presented here were derived by the computer algebra program given 
in Appendix Such a computer algebra program may be straightforwardly 
modified to model more complicated dynamical systems. 

In this work I concentrate upon a proof of the concept of applying centre 
manifold theory to constructing effective finite difference models. To that 
end we only consider an infinite domain or strictly periodic solutions in finite 
domains. Then all elements of the discretisation are identical by symme- 
try and the analysis of all elements is simultaneous. However, if physical 
boundaries to the domain of the pde are present, then those elements near 
the physical boundary will need special treatment. I do not see the need 



2 Centre manifold theory underpins the fidelity 



5 



for any new concepts, just an increase in the amount of detail. The centre 
manifold approach also sheds an interesting light upon the issue of the initial 
condition for a finite difference approximation. Earlier work on the issue of 
initial conditions in general [16, [7], |18[ hints that the initial values for the pa- 
rameters Uj(t) is not simply the field u sampled at the grid points, u(xj,0), 
but some more subtle transform. Some preliminary research suggests that a 
leading order approximation is that Uj(0) is the element average of u(x, 0), 
an approximation which usefully conserves u. Lastly, here we have only anal- 
ysed an autonomous dynamical system, the pde ([I]); forcing applied to the 
pde may be approximated using the projection obtained for initial conditions 
H [18! • Further research is needed in the above issues and in the application 



of the theory to higher order PDEs and in higher spatial dimensions. 



2 Centre manifold theory underpins the fi- 
delity 

Here I describe in detail one way to place the discretisation of PDEs within the 
purview of centre manifold theory. For definiteness I address the generalised 
Burger's equation ([!]) as an example of a broad class of PDEs. 

As introduced earlier, the discretisation is established via an equi-spaced 
grid of collocation points, Xj = jh say, for some constant spacing h. At 
midpoints Xj+1/2 = (%j + x j+i)/2 artificial boundaries are introduced with 
one extra refinement over that discussed in the Introduction: 

au+ „ (4) 

= 7 (u + - u") , (5) 

where the introduction of the near identity operator 

h 2 8 h 4 d 2 17h e d 3 

12di + 120d¥~ 201609^ + "' ' ^ 

will be explained in the next section. These boundaries divide the domain 
into a set of elements, the j'th element centred upon Xj and of width h. A 
non-zero value of the parameter 7 couples these elements together so that 
when 7 = 1 the pde is effectively restored over the whole domain. The 
generalised Burger's equation (|l|) with "internal boundary conditions" (@-|5|) 
is analysed here to give the discretisation in the interior of the domain. 
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The application of centre manifold theory is based upon a linear picture 
of the dynamics. Adjoin the dynamically trivial equation 

and consider the dynamics in the extended state space (u(x),j). This is a 
standard trick used to unfold bifurcations fl], §1.5] or to justify long- wave 



approximations [15]. Within each element u = 7 = is a fixed point. Lin- 
earized about each fixed point, that is to an error O {\\u\\ 2 + 7 2 ), the pde 
is 

du d 2 u , du 



dt dx 2 ' dx 



= 

E J±l/2 



namely the diffusion equation with essentially insulating boundary condi- 
tions. There are thus linear eigenmodes associated with each element: 



e Xnt cos[nir(x - Xj-y 2 )/h] , 2^-1/2 < x < x j+ i/ 2 
, otherwise , 



7 = 0, u oc 

for n = 0, 1, . . ., where the decay rate of each mode is 



n 2 n 2 



A„. = -^; (9) 

together with the trivial mode 7 = const, u = 0. In a domain with m 
elements, evidentally all eigenvalues are negative, —ir 2 /h 2 or less, except for 
m + 1 zero eigenvalues: 1 associated with each of the m elements and 1 from 
the trivial (0). Thus, provided the nonlinear terms in ([[]) are sufficiently 
well behaved, the existence theorem (|| p281] or |25|, p96]) guarantees that 
a m + 1 dimensional centre manifold M. exists for ([!]) with (§-0). The centre 
manifold M. is parameterized by 7 and a measure of u in each element, say 
uy. using u to denote the collection of such parameters, M. is written as 

u(x, t) — v(x; u, 7) . (10) 

In this the analysis has a very similar appearance to that of finite elements. 
The theorem also asserts that on the centre manifold the parameters Uj evolve 
deterministically 

Uj = 9j(u,l), (11) 

where Qj is the restriction of ([!]) and (^-0) to M. In this approach the 
parameters of the description of the centre manifold may be anything that 
sensibly measures the size of u in each element — I simply choose the value 
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of u at the grid points, Uj = u(xj,t). This provides the necessary amplitude 
conditions, namely that 

Uj = v\ Xj . (12) 

The above application of the theorem establishes that in principle we may 
find the dynamics (|TTD of the interacting elements of the discretisation. A 
low order approximation is given in (||]). 

The next outstanding question to answer is: how can we be sure that 
such a description of the interacting elements does actually model the dy- 
namics of the original system (|IJ) with (§-0)? In the development of inertial 
manifolds by Temam [24| and others, this question is sometimes phrased as 
one about the asymptotic completeness of the model, for example see Robin- 
son 1 22] or Constantin et al |5], Chapt.12-3]. Here, the relevance theorem of 



centre manifolds, p282] or [f25[ pl28], guarantees that all solutions of ([I]) 
and (f|-[7]) which remain in the neighbourhood of the origin in (u(x), 7) space 
are exponentially quickly attracted to a solution of the m finite difference 
equations fliTf). For practical purposes the rate of attraction is estimated by 
the leading negative eigenvalue, here —n 2 /h 2 . Centre manifold theory also 
guarantees that the stability near the origin is the same in both the model 
and the original. Thus the finite difference model will be stable if the original 
dynamics are stable. After exponentially quick transients have died out, the 
finite difference equation (pTTJ ) on the centre manifold accurately models the 
complete system ([!]) and (£| |7|)- 

The last piece of theoretical support tells us how to approximate the shape 
of the centre manifold and the evolution thereon. Approximation theorems 
such as that by Carr & Muncaster |3|, p283] assure us that upon substituting 
the ansatz (|iOH Pl|) into the original (]l|) and (§-0) and solving to some order of 
error in and 7, then M. and the evolution thereon will be approximated 
to the same order. The catch with this application is that we need to evaluate 
the approximations at 7 = 1 because it is only then that the artificial internal 
boundaries are removed. In some applications of such an artificial homotopy I 
have demonstrated good convergence in the parameter 7 [F7] . Thus although 
the order of error estimates do provide assurance, the actual error due to the 
evaluation at 7 = 1 should be also assessed otherwise. Here, as discussed in 
Section [|, we have crafted the interaction ([5]) between elements so that low 
order terms in 7 recover the exact finite difference formula for linear terms. 
Note that although centre manifold theory "guarantees" useful properties 
near the origin in 7) space, because of the need to evaluate asymptotic 

expressions at 7 = 1, I have used a weaker term elsewhere, namely "assures". 
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3 Numerical comparisons show the effective- 
ness 



We now turn to a detailed description of the centre manifold model for 
Burger's equation ([[]). 

The algebraic details of the derivation of the centre manifold model (|T0|- 
Tl~D is the task of the computer algebra program listed in the Appendix. In 
an algorithm introduced in |2T)fl , the program iterates to drive to zero the 
residuals of the governing differential equation ([3]) and its boundary condi- 
tions • A key part of the iteration is to solve for corrections v' and g' 
from the linear diffusion equation within each element 



<9V 

dx 2 



g' + R, 



dv' 
dx 



R 



± 



where R and R± denote the residuals for the current approximation. The 
initial approximation is simply that in each element u = Uj, constant, with 
iij = gj = 0. The algebraic details of this iteration are largely immaterial 
so long as the residuals iterate to zero to some order in u and 7. This is 
achieved by the listed computer algebra program, it is included to replace 
the recording of tedious details of elementary algebra. 

The shape of the centre manifold gives the field function of the 

parameters Uj and the coupling 7. To low-order in u and 7, and written in 
terms of the scaled coordinate £ = (x — Xj)/h, the solution field in the jth 
element is 

2uj + u 



u 



Uj + 7 



Uj+i - Uj-i 



auj(u j+1 - 2uj + Uj-i) ( ^f 3 



+ hy 
+ /i 2 7 

+ (3(u'j +1 - 3u j+ iUj + Au 
+ (3{u) +l - 3 Uj 



2 1l 

a Uj{Uj+i 



3 1 K 24 48 



3Un 



3-1J 
1 



24 s 48 s 



Ui 



U 



12 

2uj + Uj-.i)£ 



48 



+ o(\ 



u 



| 4 ,7 2 



(13) 



Observe the first line, when evaluated at 7 = 1, is simply the quadratic 
interpolant based upon second order accurate centred differences. This is 
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Figure 1: low order approximations (|13|) to the centre manifold field u(x) 
within each element for a fixed set of data points (circles) corresponding to 
m = 8 intervals on [0, 27r) with parameters a = 6, (3 = and 7 = 1: • ■ •, 
the normal field associated with finite differences; — includes the 
first effects of the nonlinear advection term, auu x ; — — — , includes the 
second-order effects in a. 
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normal finite differences. As displayed in Figure [I], the second line shows that 
this field should be modified because of the nonlinear advection term auu x ; 
the modification in proportion to the second difference at Xj is reasonable 
because when u has a local maximum the field must be decreased/increased 
to the left/right due to the advection to the right of the lower/higher levels 
of u respectively. Note that the same discretisation for this nonlinear term is 
obtained here whether we code it as auu x or a(u 2 ) x /2 — the centre manifold is 
independent of any valid change in the algebraic description of the dynamics. 
The third line gives the next order correction due to the nonlinear advection 
and is also displayed in Figure |I[ The fourth, fifth and sixth lines show 
how the cubic nonlinearity, (3u 3 , modifies u in the element because of its 
nonlinearly varying effect as a source or sink when the field u itself varies. It 
might appear odd to introduce the dysfunctional behaviour between elements 
shown in Figure [1|, but centre manifold theory reasonably assures us that 
it is appropriate for the nonlinear dynamics we wish to model. In short, 
the description of the centre manifold is based upon standard differencing 
formulae, but includes effects upon u within each element due to the nonlinear 
processes that act in the continuum dynamics. 

The finite difference model is given by the evolution on the centre mani- 
fold. To linear terms in 7 but to two orders higher in it is 



U; 



1 



7a 



— (u j+ i - 2 Uj + Uj-J - -^Ujiuj+i - Uj-t) - (3u 



h 2 

+ 7 



u 2 Auj + i — 2uj + Uj-i) 



a* 
12 



+ 

+ h 2 j 



U j( U j+l ~ U j-l) 



a 



720 



Uj(uj+i — 2uj + Uj-i) 



+ 



a 2 (3 
5760 

(3 2 



u) (9«J +1 - \Uu j+1 u 2 + 208^ 3 - 113«,_i«| + 9^ 



+ luj-xu) + §u)_ x u 2 - 2\u)_ x 



+^(n«ir,7 2 



(14) 



The first line recorded here, when evaluated for 7 = 1, is just the clas- 
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Figure 2: Contours of an exact solution u(x, t), , to compare with numer- 
ical approximations (15): • ■ -, the conventional approximation, errors o(a); 
— ■ — the first correction, errors o(a 2 ); — — — , the second correction, 
errors o(a 4 ). Burger's equation (HD with parameters a = 6 and (3 = is 
discretised on m = 8 elements in [0, 2ir) and drawn with contour interval 
Au = 0.2. 



sic second-order finite difference equation for the generalised Burger's equa- 
tion (|T]). The second line starts accounting systematically for the variations 
in the field u within each element and how they affect the evolution through 
the nonlinear terms. The approximation formed by the first three lines is 
that reported in the Introduction as (|3|). The fourth and further lines above, 
through a(3 and a 2 [3 effects, show that this approach also accounts for in- 
teraction between the nonlinear terms in the pde in the finite difference 
approximation. Finite difference equations derived via this approach holisti- 
cally model all the interacting dynamics of the entire pde. 

To show the effectiveness of the approach I compare the finite difference 
model QUI) to exact solutions of Burger's equation ([!]) with (3 = 0. Exact 
solutions are obtained via the Cole-Hopf transformation u = —(2/a)<f) x /<j) by 
choosing <p(x, t) satisfying the diffusion equation fa = <p xx . For example, here 
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Table 1: Comparison of the errors e, defined in ([16]), of the spatially second 
order model fll5"D at the two first truncations in the nonlinear parameter 
a for various grid subintervals m. Observe that the higher order model is 
significantly less sensitive to the nonlinearity. 





m 


= 8 


m = 


16 


m = 


32 


a 


o(a) 


o(a 2 ) 


o(a) 


o(a 2 ) 


o(a) 


o(a 2 ) 


1 


.0116 


.0095 


.0031 


.0025 


.0008 


.0006 


3 


.0356 


.0109 


.0102 


.0040 


.0026 


.0011 


6 


.0694 


.0115 


.0215 


.0059 


.0054 


.0018 


10 


.0971 


.0186 


.0318 


.0055 


.0081 


.0026 



I choose 27r-periodic functions 



14 



-e 



-7T 2 /4t 



+ ^ / exp 



(x - k2nf 
4(t + t ) 



with t = 7r /(4\/^); chosen so that the initial state u(x, 0) is a rough ap- 
proximation to sin (a;). Note that consequently \\u(x, 0)|| ~ 1. Choosing m 
intervals on [0, 27r) gives an element length h = Its jm and grid points Xj = jh 
for j = 0, . . . , m — 1. Because of the antisymmetry in u(x, t) about x = kir, 
I only display the interval [0, 7r]. With /3 = 0, there are three different ap- 
proximations from (|14D , with 7 = 1, depending upon where the expansion is 
truncated in ||w|| (or equivalently in a): 



1 a 



a n/ \ 



ct 4 

-h 2 —u 4 j {u j+l -2u j + u j - 1 ). (15) 

Just the first line forms a model with o(a) errors (a normal finite difference 
approximation), the first two lines form a model with o(a 2 ) errors, and all 
shown terms form the model with o(a 4 ) errors. The solutions of these models 
over < t < 1 with m = 8 and nonlinearity parameter a = 6 are shown in 
Figure [| Observe that the leading approximation (dotted) is significantly 
in error whereas the next two refinements (dot-dashed and dashed) are rea- 
sonably accurate. Such accuracy is remarkable considering the high level of 
nonlinearity, a = 6, and the few points in the discretisation, m — 8. 

For an assessment of errors over a range of parameters I record in Table |I| 
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the error 

e = m f x f^EI^(*)lj > ( 16 ) 

where (t) is the difference between an approximation and the exact solution 
at the jth grid point. Observe the usual properties that the errors decrease 
with increasing spatial resolution m, and that the errors increase with in- 
creasing nonlinearity a. Our interest lies more in the comparison between 
the o(a) and o(a 2 ) models (here the o(a 4 ) model is virtually indistinguish- 
able from the o(a 2 ) model). For small nonlinearity, a = 1, there is very little 
difference; presumably the errors are dominated by the second-order approx- 
imation of spatial derivatives. As the nonlinearity a increases, the error of 
the o(a) approximation increases roughly in proportion to a, but the error 
of the o(a 2 ) approximation grows much more slowly. We conclude, as centre 
manifold theory assures us, that the nonlinearly corrected o(a 2 ) model more 
accurately captures the nonlinear dynamics of Burger's equation. 



4 Higher order approximations are more ac- 
curate 

So far I have reported results only to first order in the coupling coefficient 7. 
Retaining higher orders in 7 gives higher order difference rules as the coupling 
between adjacent elements is always ameliorated by a power of 7. However, 
it is only with the specific coupling of @ that the 2p + 1 width stencil, 
obtained by retaining terms in / y p , attains 2pth-order accuracy in the linear 
terms. Retaining j 2 terms, for example, then gives a fourth-order model 
which also performs remarkably well, especially for larger nonlinearities. 

The particular form of the artificial internal boundary condition ([5]) con- 
trols the flow of information between elements of the discretisation. When 
7 = 1 the condition reduces to requiring the desired continuity: u~ = u + 
from the rightmost term and u~ = w+ from (|]) where, as before, the su- 
perscripts ± denote evaluation at the internal boundary of expressions from 
the right /left-hand elements respectively. But we have enormous freedom in 
choosing the (1 — 7)-term in ([5]). Our first requirement is that when 7 = 
it insulates the elements from each other so that the values of Uj in each 
element are independent in the centre manifold analysis. This is achieved is 
conjunction with (f|) by only involving u~+u^ and its time derivatives. This 
is all that is necessary to apply centre manifold theory. 

However, computer algebra experiments show that generally we have to 
sum to high-order in 7 to recover standard finite difference formulae for the 
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linear diffusion term u xx . This is impractical because the width of the finite 
difference formula also grows with the order of 7. But we have freedom to 
find an interaction where the expansions of the linear terms truncate in 7 
and thus recover the standard formulae at low-orders in 7. Computer algebra 
experiments show that the expansion (||) for A ensures the linear diffusion 
term is modelled with errors O (h 2p ) by a 2p + 1 width finite difference stencil 
when the centre manifold is constructed with errors O (7 P+1 ), for p = 1, . . . ,4. 
The expression for A in (|6]) is recognised to be the "asymptotic expansion" 
of the operator 



It is not apparent why this particular operator is so desirable. However, 
an answer may lie via the following observation. In the linear dynamics, 
Ut = u xx and so effectively d t = d x . Denoting the fundamental spatial 
differential operator T> = hd x , observe that in effect 



where 5 denotes the centred difference operator and fi the centred averaging 
operator: 



Using this in (|), "cancelling" V and "multiplying" by /i leads to the internal 
boundary condition being equivalent to 



This has a pleasing symmetry as (u + + u~)/2 on the left is like the averaging 
operator fi on the right, and the (u + — u~) on the right is like the differencing 
operator 5 on the left. However, at this stage I also do not know why this 
should be desirable. I leave this aspect for further research. 

Approximations to fourth-order accuracy in space will be explored to 
verify their accuracy. Fourth-order approximation formulae are obtained by 
truncating the analysis to errors O (7 s ). This enables the direct interaction 
between neighbouring elements through the O (7) terms and the indirect 
interaction with the next nearest neighbours through O ( r y 2 ) terms. Running 
the computer algebra program listed in Appendix [A] to higher order in 7 
gives the following for the evolution of the amplitudes Uj (using the centred 




2 

A= — tanh(X>/2) 



8 




(1 — 7)<5-(w + + u ) = 7/i(-u + — u ) 



(17) 
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Table 2: Comparison of the errors e, defined in (H^), of the spatially fourth 
order model fllTf ) for the first and third truncations in the nonlinear parameter 
a for various grid subintervals m. Observe that the higher order model is 
significantly less sensitive to the nonlinearity. 





m 


= 8 


m = 


16 


a 


o(a) 


o(a 3 ) 


o(a) 


o(a 3 ) 


1 


.0013 


.0016 


.0003 


.0003 


3 


.0070 


.0090 


.0009 


.0021 


6 


.0252 


.0131 


.0020 


.0052 


10 


.0464 


.0132 


.0066 


.0066 


20 


.0659 


.0171 


.0179 


.0056 


30 


.0683 


.0196 


.0250 


.0061 
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The first line are terms from the model discussed in the previous section, 
but here shown to lower order in \\u\\. The second line gives the appropriate 
correction for the linear diffusion term u xx to make it fourth-order accurate in 
space. The third line gives a fourth-order correction for the advection terms 
auu x . The fourth and further lines give spatially higher-order corrections to 
the nonlinearly higher-order terms. 

These fourth-order in space terms ensure a higher accuracy. Setting the 
artificial parameter 7 = 1 and truncating at various orders in \\u\\, equiva- 
lently in a because I set ft = as before, we find that to see any appreciable 
difference between the nonlinearly higher-order models we have to increase 
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Figure 3: Contours of exact solution u(x, t), , to compare with fourth- 
order numerical approximations (|i8|): • • • •, the conventional approximation, 
errors o(a); — • — the first correction, o(a 2 ); and — — — includes the 
second correction, o(a 3 ). Burger's equation ([[]) with parameters a = 10 and 
(3 = is discretised on m = 8 intervals and drawn with contour interval 
Am = 0.2. 
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the nonlinear parameter to a = 10 or more (from the 6 used in the previous 
section). See in Figure ^| that the nonlinearly low-order approximation is a 
little in error, but that the higher-order models are generally better. They 
are reasonably accurate despite the large nonlinear ity, a = 10, even though 
there is only m = 8 intervals in one period. A more comprehensive summary 
of the numerical errors in given in Table |2] which compares the overall error 
e for m = 8 and m = 16 elements over one period for the nonlinearity pa- 
rameter ranging over 1 < a < 30. Once again observe that the conventional 
o(a) model has errors roughly proportional to a, whereas the o(a 3 ) model is 
much less sensitive. As the model is fourth-order in space the errors here are 
an order of magnitude smaller than those in Table p] of the previous section. 
This confirms the effectiveness of this approach to developing finite difference 
models. 

However, because of the subtleties in the analysis of the nonlinear terms, 
the equivalent partial differential equation of a derived finite difference model 
will not reduce to the original pde to the expected order in h. For example, 
here the equivalent pde of the model (0), with 7 = 1, is 

du du d 2 u _ , h 2 / n , ,\ ,./, 4 \ 

dt + aU d~x = ~ PU + 12 \- 2au * u ™ + a uu *) + ° \ h ) ' ( 19 ) 

which apparently differs from the generalised Burger's equation by 0(h 2 ). 
But note that this equivalent differential equation is obtained as h — > 
keeping ||w|| fixed, whereas the centre manifold model (|1^) is derived for 
fixed h as ^0 but taking full account of nonlinear dynamics within the 
domain. The numerical results presented here support my contention that 
this centre manifold approach better models the pde for finite h and 



5 Conclusion 

Using centre manifold theory is a powerful new approach to deriving finite 
difference models of dynamical systems. Many details need to researched 
for a general application of the theory. However, this particular example 
application to a generalised Burger's equation ([!]) shows many promising 
features. 

By introducing artificial internal boundaries we apply centre manifold 
theory (§0). The internal boundaries divide the domain of interest into sub- 
domains that look very much like finite elements. The theory guarantees 

x It is conceivable that a more carefully crafted interaction between the elements, 
modifying (||), may achieve O (h 2p ) consistency between the finite difference model 
and all terms of the original pde when constructing such a centre manifold model to 
0(~/P+\\\up). 
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essential properties required for a low- dimensional model, though the need to 
evaluate the asymptotic expressions for 7 = 1 weakens this assurance. First, 
a model exists parameterized in any reasonable way we desire. Second, the 
model is approached exponentially quickly by the original dynamics, in a time 
of O (h 2 ) for Burger's equation, but possibly much quicker for higher-order 
PDEs. The same theorem also assures us that the numerical model shares the 
same stability as the original pde dynamics. Lastly, an approximation to the 
model may be found to any order in the nonlinearity, using computer algebra 
for example. The resultant finite difference models are independent of valid 
rearrangements of the governing pde. In effect this technique analyses the 
actual dynamics of the pde as a whole — this is a holistic approach. For 
example, I am not aware of any other discretisation method that develops 
cross terms between the nonlinearities, yet the presence of af3 and a 2 /3 terms 
in flH|) suggests that they are essential for good finite difference models. 
Specific numerical simulations in and §(| show that the finite difference 
models derived for Burger's equation are indeed accurate. 

Further research is needed to incorporate physical boundary conditions 
into the model. The extra analysis should be straightforward, but the level 
of detail would increase significantly as near the physical boundary we would 
lose the translational symmetry between the elements. Then the intriguing 
issue of initial conditions needs further work. Following |H| [7], [TS[ we note 
that to make accurate forecasts with the numerical models we need to provide 
initial values u® which are not the initial field uq(x) sampled at the collocation 
points Xj. Instead preliminary work shows that should be approximately 
the average of Uq(x) over each element. Non-autonomous forcing of the pde 
will need projecting onto the finite difference grid in a similar manner. The 
analysis of other PDEs in possibly higher-spatial dimensions would appear to 
hinge upon being able to solve the linear part of the pde over the elements 
subject to forcing from the coupling across the artificially introduced internal 
boundaries. This is potentially quite complicated, but the results of the 
simple problems here suggest that the analysis could be very worthwhile. 



A Computer algebra develops the approxi- 
mations 

Just one of the virtues of this centre manifold approach to modelling is that 
it is systematic. This enables relatively straightforward computer programs 
to be written to find the centre manifold and the evolution thereon p(| e.g.]. 



I believe computer algebra will be increasingly used to support research by 
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performing extensive routine algebraic manipulations. To ensure correctness 
and to provide a basis for further work it seems reasonable to include the 
computer algebra code. 

For this problem the iterative algorithm is implemented by a computer 
algebra program written in reduce | Although there are many details in 
the program, the correctness of the results are only determined by driving to 
zero (line 59) the residual of the governing differential equation, evaluated on 
line 48, to the error specified on line 46 and with the residual of the internal 
boundary condition computed on lines 50-53. The other details only affect 
the rate of convergence to the ultimate answer. 

1 COMMENT Try holistic finite differences on a generalised Burgers 

2 equation. Consider ODEs of the form u_t - u_xx = f(u,u_x) for 

3 strictly nonlinear f where for Burgers equation f \propto -u.u_x . 
4 

5 Discretise the problem by mathematically introducing insulating 

6 boundaries at grid midpoints, and coupling them together via a 

7 parameter \gamma so that the original C~2 holds at the midpoints when 

8 \gamma=l. Treat \gamma as a perturbation parameter. 
9 

10 The unknown u(x,t) field is then a function of u_j (t)=u(x_j ,t) , where 

11 each u_j evolve according to du_j/dt=g_j . 
12 

13 Tony Roberts, 12 Sept 1998; 
14 

15 % improve printing 

16 on div; off allfac; on revpri; factor gam,usz,h; 
17 

18 7o make function of xi=(x-x_j)/h 

19 depend xi,x; 

20 let df (xi,x)=>l/h; 
21 

22 7 parameterise with evolving u(j) 

23 operator u; depend u,t; 

24 let df (u(~k) ,t)=>sub(j=k,g) ; 
25 

26 7, intx computes \int_{-l/2}~{+l/2} . . . dxi 

27 operator intx; linear intx; 

28 let { intx(xi~~p,xi) => (l+(-l)~p)/(p+l)/2~ (p+1) 

29 , intx(xi,xi) => 

30 , intx(l,xi) => 1 }; 
31 

32 7. solves v_xixi=RHS s.t. v(0)=0 and v_xi(+l/2)+v_xi(-l/2)=0 

33 operator solv; linear solv; 

2 At the time of writing, information about reduce was available from Anthony 
C. Hearn, RAND, Santa Monica, CA 90407-2138, USA. E-mail: reduce@rand.org 
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34 let { solv(xi~~p,xi) => (xi~ (p+2) /(p+2) -(1- (-1) ~p) *xi/2~ (p+2) )/(p+l) 

35 , solv(xi.xi) => (xi~3-3*xi/4)/6 

36 , solv(l,xi) => (xi~2)/2 }; 
37 

38 7 linear solution in jth interval 

39 % usz is used to measure the order of terms in u_j 

40 v:=usz*u(j); 

41 g:=0; 
42 

43 7, iteration 

44 7. the BC is crafted s.t. gam~p kept gives (2p)th order in space 

45 7. by (l-gam)*delta(vp+v)h/2-gam*mu(vp-v)=0 recognising u_t=u_xx 

46 let { gam~3=>0, usz~4=>0 } ; 

47 repeat begin 

48 deq:=df (v,t)-df (v,x,2)+(a*v*df (v,x)+b*v~3) ; 

49 vp:=sub(j=j+l,v) ; 

50 dv:=sub(xi=-l/2,vp)-sub(xi=l/2,v) ; 

51 svd:=(h/2)*(sub(xi=-l/2,df (vp,x))+sub(xi=l/2,df (v,x))) ; 

52 bcp:=(l-gam)*( svd -df (svd,t)*h~2/12 +df (svd,t,2)*h~4/120 

53 -df (svd,t,3)*h~6*17/20160 ) -gam*( dv ); 

54 bcm:=sub(j=j-l,bcp) ; 

55 gd:=(-bcp+bcm)/h~2-intx(deq,xi) ; 

56 g:=g+gd/usz; 

57 v : =v+h~2*solv (deq+gd , xi) -xi* (bcm+bcp) /2 ; 

58 showtime; 

59 end until (deq=0) and(bcp=0) ; 

60 bcp:=sub(xi=-l/2,df (vp,x))-sub(xi=l/2,df (v,x)) ; 7. check other IBC 

61 deq:=sub(xi=0,v)-usz*u(j) ; % check amplitude 
62 

63 end; 
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